$ontext
This code is for the following paper:
    Cai, Yongyang, William Brock, and Anastasios Xepapadeas (2022).
    Climate Change Impact on Economic Growth: Regional Climate Policy under Cooperation and Noncooperation.
    JAERE, forthcoming

Author: Yongyang Cai, The Ohio State University
Date: August 11, 2022
$offtext

parameters
    a1N
    a2N
    a3N
;    

set ii /North, Tropic, South/;
set r(ii) /North, Tropic/;

* years in CMIP5 data
set tC / 2006*2100 /;

* years in SSP data
set t / 2015*2515 /;
set tfirst(t);
tfirst(t)$(ord(t)=1) = YES;
alias(s,t);

set gd /ga0, dela/;

$offlisting

* temperature anomaly for RCP from CMIP5 in 2010-2100

Table TA(tC,ii) 
$ondelim
$include RegionalSurfaceTempAnomaly_RCP26.csv
$offdelim
;

Table Pop(t,r) 
$ondelim
$include RegionalPop_SSP1_Annual.csv
$offdelim
;

Table paraTFP(gd,r) 
$ondelim
$include paramTFP_noDam.csv
$offdelim
;

$onlisting

parameters
    gama     Capital elasticity in production function        /.300    /
    dk       Depreciation rate on capital (per year)          /.100    /

    elasmu   Elasticity of marginal utility of consumption     /1.45 /
    prstp    Initial rate of social time preference per year   /.015  /

    ga0(r)      Initial growth rate for TFP per year             
    dela(r)     Decline rate of TFP per year
    A0(r)       initial TFP
    ga(t,r)
    A_noDam(t,r)
    A(t,r)      TFP    

	k0(r)
	GDP0(r)
	betas(t)   compound discount factor
	
	TA0(r)
    TATM(t,ii)	
	TAdev(t,r)
;

ga0(r) = paraTFP('ga0',r);
dela(r) = paraTFP('dela',r);

betas(t) = 1/(1+prstp)**(ord(t)-1);

* In Burke et al. (2018), GDP in 2010 is the multiplication of population in 2010 and
* the average of GDP per capita from 1980 to 2010 (after dropping those items with NA).
* However, I need to set the initial year to be 2015, so I use world bank data
* This will only affect the value of initial TFP, as we are just matching the growth rates of GDP

* Use world bank data (2015 GDP) 
GDP0('North') = 54;
GDP0('Tropic') = 19;

* RICE data in 2015
k0('North') = 100;
k0('Tropic') = 53;

TA0(r) = sum(tC$(tC.val=2015), TA(tC,r));

TATM(t,r)$(t.val<=2100) = sum(tC$(tC.val=t.val),TA(tC,r));
TATM(t,r)$(t.val>2100) = TATM('2100',r);

TAdev(t,r) = TATM(t,r) - TA0(r);

parameters
    a1(r)
    a2(r)
    a3(r)
;    

a1('North') = %a1N%;
a2('North') = %a2N%;
a3('North') = %a3N%;
a1('Tropic') = 0;
a2('Tropic') = 0;
a3('Tropic') = 0;

ga(t,r)=ga0(r)*exp(-dela(r)*(ord(t)-1));

A0(r) = sum(tfirst, GDP0(r) / (Pop(tfirst,r)**(1-GAMA)*k0(r)**GAMA));
A_noDam(tfirst,r) = A0(r);
loop((t,r)$(ord(t)<card(t)),
  A_noDam(t+1,r) = A_noDam(t,r)/(1-ga(t,r));
);

A(t,r) = A_noDam(t,r)/( 1+sum(s$(ord(s)<=ord(t) and ord(s)>=ord(t)-10),(a1(r)*TAdev(s,r)+a2(r)*sqr(TAdev(s,r)))*a3(r)**(ord(t)-ord(s))) );
display A;

positive variables
 K(t)            Capital stock trillions US dollars
 Y(t)
 c(t)
;
c.lo(t) = 0.01;
K.lo(t) = 0.01;

VARIABLES
 obj
;


EQUATIONS
 KK(t)
 GrossOutput(t)

 objfun
;

KK(t)$(ord(t)<card(t))..       K(t+1) =e= (1-DK)*K(t) + Y(t) - c(t)*Pop(t,'North');
GrossOutput(t)$(ord(t)<card(t))..	Y(t) =e= A(t,'North')*Pop(t,'North')**(1-GAMA)*K(t)**GAMA;

Objfun..            obj =e= SUM( t$(ord(t)<card(t)), betas(t)*c(t)**(1-elasmu)/(1-elasmu)*Pop(t,'North') );

K.fx(tfirst) = k0('North');

K.l(tfirst) = k0('North');
loop(t$(ord(t)<card(t)),
  Y.l(t) = A(t,'North')*Pop(t,'North')**(1-GAMA)*K.l(t)**GAMA;
  c.l(t) = 0.7*Y.l(t)/Pop(t,'North');
  K.l(t+1) = (1-DK)*K.l(t) + Y.l(t) - c.l(t)*Pop(t,'North');
);  

**********************************

** Solution options
option iterlim = 99900;
option reslim = 99999;
option solprint = off;
option limrow = 0;
option limcol = 0;
option nlp = conopt;


model Growth / all /;

solve Growth maximizing obj using nlp ;

ABORT$(Growth.MODELSTAT>2 or Growth.SOLVESTAT>1) "Model not normally completed";

parameter capY(t);
capY(t) = Y.l(t) / Pop(t,'North');

file results /output_Dam_lag10_North.csv/;
results.nd = 6 ;
results.nw = 0 ;
results.pw=20000;
results.pc=5;

put results;

loop(t$(t.val<2100),
  put capY(t):14:6; 
  put /;
);



